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Abstract 

During the past decades the study of strongly interacting fluids experienced a tremendous 
progress. In the relativistic heavy ion accelerators, specially the RHIC and LHC colliders, it 
became possible to study not only fluids made of hadronic matter but also fluids of quarks and 
gluons. Part of the physics program of these machines is the observation of waves in this strongly 
interacting medium. From the theoretical point of view, these waves are often treated with li- 
nearized hydrodynamics. In this text we review the attempts to go beyond linearization. We show 
how to use the Reductive Perturbation Method to expand the equations of (ideal and viscous) rela- 
tivistic hydrodynamics to obtain nonlinear wave equations. These nonlinear wave equations govern 
the evolution of energy density perturbations (in hot quark gluon plasma) or baryon density per- 
turbations (in cold quark gluon plasma and nuclear matter). Different nonlinear wave equations, 
such as the breaking wave, Korteweg-de Vries and Burgers equations, are obtained from different 
equations of state (EOS). In nuclear matter, the Walecka EOS may lead to a KdV equation. We 
explore equations of state such as those extracted from the MIT Bag Model and from QCD in the 
mean field theory approach. Some of these equations are integrable and have analytical solitonic 
solutions. We derive these equations also in spherical and cylindrical coordinates. We extend the 
analysis to two and three dimensions to obtain the Kadomtsev-Petviashvili (KP) equation, which is 
the generalization of the KdV. The KP is also integrable and presents analytical solitonic solutions. 
In viscous relativistic hydrodynamics we have second order patial derivatives which physically rep- 
resent dissipation terms. We present numerical solutions and their corresponding algorithms for 
the cases where the equations are not integrable. 



I. INTRODUCTION 



The elementary particles and their interactions are well described by the Standard Model 
(SM), which is an extremely successful theory [T]. In this theory matter is composed by 
quarks and leptons and their interactions are due to the exchange of gauge bosons. The 
sector of the SM which describes the strong interactions at the fundamental level is called 
Quantum Chromodynamics, or QCD [2]. According to QCD there are six types, or flavors, 
of quarks: up, down, strange, charm, bottom and top, and they interact exchanging gluons. 
Quarks and gluons have a special charge called color, responsible for the strong interaction. 
They do not exist as individual particles but, due to the property of color confinement, quarks 
and gluons form clusters called hadrons, which can be grouped in baryons and mesons. The 
former are made of three quarks, as the proton, and the latter are made of a quark and an 
antiquark, as the pion, for example. The quarks carry a fraction of the elementary electric 
charge and a fraction of the baryon number. The electric and color charges and the baryon 
number are conserved quantities in QCD. 

The hadrons also interact strongly and their interactions have been traditionally described 
by a quantum field theory called Quantum Hadrodynamics, or QHD, in its different versions 
[2111]. According to QHD, in nuclear matter neutrons and protons interact exchanging scalar 
and vector mesons. This potentially complicated theory becomes quite simple in the mean 
field approximation, in which the meson fields are treated as classical fields. The different 
versions of QHD in the mean field approximation are called relativistic mean field models 
(RMF)j3j Hj. More recently, nuclear matter has been studied with Effective Field Theories, 
which incorporate the fundamental symmetries of QCD in hadron physics. 

In the phase diagram of QCD, we can observe that under extreme conditions of very 
large temperatures and/or very large densities, the normal hadronic matter undergoes a 
phase transition to a deconfined phase, a new state of matter called the quark gluon plasma 
(QGP) [5]. Together with deconfinement, a second phase transition takes place: the chiral 
phase transition, during which chiral symmetry is restored and the light quarks (up and 
down) become massless. The hot QGP is produced in relativistic heavy ion collisions in 
the Relativistic Heavy Ion Collider (RHIC) at the Brookhaven National Laboratory (BNL) 
[6j [7] and even more in the Large Hadron Collider (LHC) at CERN. The cold QGP may 
exist in the core of compact stars |8]. 
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According to our present understanding of the RHIC measurements, the QGP behaves as 
an almost perfect fluid and its space-time evolution can be very well described by relativistic 
hydrodynamics [5HTT]. The discovery of this new fluid motivated inumerous theoretical 
works addressing viscosity in relativistic hydrodynamics [T2HH] . At the same time, more 
sophisticated measurements made possible to study the propagation of perturbations in the 
QGP. We may, for example, study the effect of a fast quark traversing the hot QGP medium. 
As it moves supersonically throughout the fluid, it may generate waves of energy density or 
baryon density [15J. It may be even possible that these waves may pile up and form Mach 
cones, which would affect the angular distribution of the produced particles, fluid fragments 
which are experimentally observed. 

The study of waves in the quark-gluon fluid has been mostly performed with the assump- 
tion that the amplitude of the perturbations is small enough to justify the linearization of 
the Euler and continuity equations [1 1] . The analysis of perturbations with the linearized 
relativistic hydrodynamics leads to the standard second order wave equations and their 
traveling wave solutions, such as acoustic waves in the QGP. While linearization is justified 
in many cases, in others it should be replaced by another technique to treat perturbations 
keeping the nonlinearities of the theory. Since long ago there is a technique which preserves 
nonlinearities in the derivation of the differential equations which govern the evolution of 
perturbations. This is the reductive perturbation method (RPM) [ToT - [T9] . 

Nonlinearities may lead, as they do in other domains of physics, to new and interesting 
phenomena. In a pioneering work |20j, with the use of nonrelativistic ideal hydrodynamics 
combined with the RPM and with an appropriate equation of state of cold nuclear matter, it 
was shown that it is possible to derive a Korteweg- de Vries (KdV) equation for the baryon 
density, which has analytic solitonic solutions. This suggests that a pulse in baryon density 
(the KdV soliton) can propagate without dissipation through the nuclear medium. If this 
occurs this pulse might be responsible for an interesting phenomena: the apparent nuclear 
transparency in proton nucleus collisions at low energies. The incoming proton would be 
absorbed by the fluid (target nucleus) and turned into a density pulse. In the fluid the pulse 
satisfies the KdV equation and it is able to traverse the target nucleus emerging on the other 
side. This can be called "nuclear transparency" and it is illustrated in the Fig[T} 

Perturbations in fluids with different equations of state (EOS) generate different nonlinear 
wave equations: the breaking wave equation, KdV, Burgers... etc. Among these equations 
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FIG. 1: KdV soliton in the baryon density simulating "nuclear transparency". 



we find the Kadomtsev-Petviashvili (KP) equation [21], which is a nonlinear wave equation 
in three spatial and one temporal coordinate. It is the generalization of the KdV equation 
to higher dimensions. This equation has been found with the application of the reductive 
perturbation method [IB] to several different problems such as the propagation of solitons 
in multicomponent plasmas, dust acoustic waves in hot dust plasmas and dense electro- 
positron-ion plasma [22H34"] . Previous studies on nonlinear waves in cold and warm nuclear 
matter can be found in [20~1 l35H4*0"] . Works on nonlinear waves in cold QGP in the mean 
field approach were published in [H] and their extension to three dimensions was published 
in 02]. 

In this text we review the applications of the RPM [To r tl2"2"H55] to relativistic fluid dynamics 
[9j [TO]. In the next section we review the basic equations of relativistic hydrodynamics. In 
section III we present the RPM and give the coordinate transformations to be used in the 
subsequent sections. In section IV we introduce the equations of state of hadronic matter and 
of the quark-gluon plasma, giving special attention to the ingredients which may generate 
solitons. In section V we list the wave equations which follow from the application of the 
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RPM to hydrodynamics. In section VI we present the analytical solutions of some of these 
differential equations and in section VII we show the time evolution of several nonlinear 
waves obtained by numerical integration. In section VIII we make some final remarks. 

II. RELATIVISTIC HYDRODYNAMICS 
A. Definitions and basic equations 

As we have mentioned before, the hot and dense medium created in heavy ion collisions 
at RHIC behaves approximately as a perfect fluid and ideal hydrodynamics can applied to 
describe its space-time evolution. Moreover the study of perturbations in the fluid, such as 
the waves created by fast partons, can be studied in the context of hydrodynamics as well. 
In this section we present the formalism of relativistic hydrodynamics. Pedagogical texts on 
this subject can be found in [HI HO]- Here we give special attention to the variables which 
are relevant for strongly interacting fluids, such as the baryon density. We write the final 
equations in a more extended form, which allows the reader to make a direct comparison 
with the corresponding non-relativistic versions of these equations. In what follows we use 
c = 1, ti — 1 and the Boltzmann constant is taken to be one, i.e., ks = 1. All the relativistic 
equations are written in terms of 4- vectors and the metric tensor is given by g^, with 
goo = —gn = —g22 = —g33 = 1 and g^ — if // ^ v. From [HI QUI EE] the fundamental 
equations of a relativistic ideal fluid are: 

De + (e + p)9X = ° W 

and 

(e + p)Du a - V> = (2) 

where: 

D = «"9„ and V a = A^ (3) 

and the projection operator on the orthogonal direction to the fluid velocity is given by 
the tensor: 

A"" = f v - u^u" (4) 
which has the properties A^u^ = A^ u u u = and A^A" = A^ a . 
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The velocity 4- vector of the fluid element is given by [9j [TOj, [12]: = (7,7??) where 7 
is the Lorentz factor 7 = (1 — t> 2 ) -1 / 2 and thus u^u^ = 1 and u^d v u^ = (l/2)d u (u^u IM ) = 
(1/2)^(1) = 0. 

In an ideal fluid all dissipative (viscous) effects are neglected and in order to introduce 
the effects of viscosity, it is necessary to add the viscous stress tensor IP 1 * to the energy- 
momentum tensor. In a simple way we may write for a viscous fluid: 

T*" = Tffi + if" (5) 

where 

Tjfi = euV - pA» v (6) 

is the ideal relativistic fluid energy-momentum tensor [HI [TO] . We also consider for simplicity 
[12] a system without conserved charges (or at zero chemical potential) and so the total 
momentum density is due to the flow of energy density ii M T A " / = eu v and hence we must 
assume that: 

u^W v = (7) 

We take the appropriate projections of the conservation equations of the energy momentum 
tensor: the parallel {u v d^T^ u ) and perpendicular (A^d^T^) to the fluid velocity The results 
are: 

u v d^ v = De + (s + p)dX + u v d^W v = (8) 

and 

A^T"" = (e + p)Du a - V> + A^IF" = (9) 
Using the symmetrization notation: 

A {ll B v) = l -(A^B v + A v B^) (10) 



we are able to rewrite the u u d l Jl tMl/ term in Q as u l/ d fJ Tl tlL ' = d^Uyli^) — Yi^d^Uy). We 
also use the identity 

d lt = u /t D + V lt (11) 

and the choice of frame w^IP 1 ' = 0. The fundamental equations of relativistic viscous fluid 
dynamics are finally given by: 

De + (e + p)dy - n^V(^) = (12) 



and: 



(e + p)Du a - V> + A^IP" = 



(13) 



B. The viscous stress tensor 

The viscous stress tensor IP" has not been specified yet. We shall derive it with the help 
of the second law of thermodynamics, which states that entropy density must always increase 
locally. Moreover we shall consider a system in thermodynamic equilibrium (dp = 0) with 
zero chemical potential. Based on these statements, the thermodynamic relations are given 
by: 

£ + p = Ts and Tds = de (14) 

The second law of thermodynamics is considered in the covariant form: 

> (15) 

where the 4-current s M is given by: 

s" = su^ (16) 



Inserting relations (14) into the second law (15) and assuming that Dp = we find: 



= Ds + sdX = ~ + ^f^X ( 17 ) 



Using now equation (12) we obtain: 



= in^v^j > o (is) 

The viscous tensor may be decomposed into a traceless part tt' 1u (tt 1 ^ = 0) and a non-vanishing 
trace part, so that: 

IF" = tt^ + A^n (19) 
In terms of the notation introduced in [T2] : 

= 2V^u v) - ^A^V a w" (20) 

it can be shown that: 

u^ v = (21) 

uW {li u v) = (22) 



<rv ( ^> = o 



and 



A^A M „ = 3 



Inserting (19) to (24) into (18) we find 



= -L^v^ + ^nv Q w a > o 

and this inequality (a positive sum of squares) is satisfied if: 

tt^ = r]V^u u) , II = (V a u a , r] > and ( > 



(23) 
(24) 

(25) 
(26) 



where rj is the shear viscosity coefficient and ( is the bulk viscosity coefficient. The final 



expression for the viscous tensor is given by the substitution of (26) in (19): 

IF" = r)V {ti u v) + (A^ u V a u a 



(27) 



C. The relativistic Navier-Stokes equation 



The system of equations (12), (13) and (27) is called relativistic Navier-Stokes equation. 



The temporal component (a = 0) of (13) multiplied by v l is given by: 



(e + p)(Di)v* - v l V°p + i/A°«9 M IF" = 



and the spatial component (a = i) of (13): 

(e+p^DyW = -(e + pfriDv*) + V l p - At^IP" 



Inserting (29) into (28) we find: 



(e + ph(Dv l ) + ( V *V U - V> - (VA" - A^n 
which can be rewritten with the help of (|3| and Q as: 



T/1U 



(28) 



(29) 



(30) 



d 



A dp 



""' - (31) 

The ideal relativistic fluid limit IF" = (rj = ( = 0) in the last equation provides the 



relativistic version of the Euler equation 0QII1II2]- For future use we rewrite (31) in detail. 
To do so, we recall (]3]), Q and (27) to obtain the expression for c^IF" in (31): 



<9 M IF" = V \ d^u u + d^v? - d„ 



7^ + i?-V )(//'' 



#7 
dt 



+ V • (~/v) 



di 



(32) 



and then insert this last result into (31). We find after some algebra (13 



£ + p)rf ( ^- + v ■ v + v 



dt 



dp 
di 



r]v\ 8^0^ + d. 



du^ 
~dt 



d„ 



d 



d^ 
dt 



+ V • (tu) 



^7 
<9t 



+ V • (tu) 



+), 9/( T iJ)-W-9„ 



4 i ''^) ( " r *"" ! 



0' 



C-3^)V 



^7 



+ V ■ (7#) 



+ V • (7u) 



(33) 



which is the relativistic version of the Navier-Stokes equation. The perfect fluid, described 



by i] = C = in (33) gives the relativistic version of Euler equation [HI ITU] : 

dv 



dt 



+ (v-V)v 



1 / - _^dp 



(34) 



D. Causality and the relativistic Navier-Stokes equation 



The relativistic Navier-Stokes equation (33) does not constitute a causal theory 



This 



fact can be understood when small perturbations are considered in a system in equilibrium 
with energy density Eo, pressure po and at rest (v = 0). 
Such perturbations can be described by: 



e = Eo + <fe(t, x) 



p = p + 8p(t, x) and u 1 " = (1, 0) + 5u^{t, x) (35) 



where " 5 " denotes small deviation from equilibrium. For simplicity we consider perturba- 
tions which depend only on the x space coordinate. For the particular direction a = y we 



insert (35) into (13) to find 



{e + 5e + Po + 5p)D(5u y ) - V y ( Po + 8p) + A^IP" = 



Considering only the dependence on the x coordinate it becomes: 



d 



d 



£o+Po)^ + -n^ + o(^) = o 



dt 



dx 



(36) 
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Analogously we insert (35) into the viscous tensor (19) for the particular direction v = y 



and considering only the x coordinate dependence: 







ox 



Substituting (37) in (36) we find: 



(e 0+P0 )^5u*- V -^5uy = O(5 2 ) 



(37) 



(38) 



which provides, after performing the linearization approximation (i.e., neglecting 0(5 2 )) the 
diffusion-type evolution equation: 



d n d 2 

— Su y - — 5u y = 

dt (e + p ) dx 2 



(39) 



for the perturbation 5u y = 5u v (t,x). 

The acausality property of the Navier-Stokes equation can be studied from the Laplace- 
Fourier wave ansatz for 5u y : 



5u y (t,x) = f( w ,k) 



Akx—wt 



(40) 



where f( w ,k) are the coefficient for a given frequency w and a given wave number k. When 



inserting (40) in (143) we obtain the dispersion-relation of the diffusion equation: 

V 



w 



(eo + Po) 

which gives the speed of diffusion of a mode with wavenumber k : 



dk (e + po) 



The causality violation occurs when 



lim V di f f (k) -)> oo 

fc— >-oo 



(41) 



(42) 



(43) 



which exceeds the speed of light, i.e , the speed of diffusion grows without bound for suffi- 
ciently large wavenumber due its linear dependence on the wavenumber. 

A possible way to regulate the viscous fluid theory is given by the "Maxwell-Cattaneo 
law". In this approach a new transport coefficient called relaxation time (r n ) is added to 



equation (37) as follows 



7V— IP* + IT* = -ri— 5u y + 0{5 2 ) 
at ox 



(44) 
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The derivative of the above equation with respect to x provides: 

d 2 



d d d 
t^— — Tl xy + — U xy 
w dt dx dx 



-n 



dx 2 



5u y + 0(S 2 



(45) 



which becomes, after substituting J^n xy from (36) and performing the linearization approx- 
imation, the following evolution equation: 



d 2 







6u y 



V 



d 2 



5u y = 



(46) 



dt " (e + Po) dx 2 

where the term with second derivative with respect to time provides the following dispersion 



relation when (40) is used 
1 



w 



± 



1] k 2 



2t w \ At^ 2 (e +p )T n 



\w\ 



WW 



k' 2 



— ±ik 11 
2r w V ( £ o + Po) r * 

V 



1 



Ar 2 k 2 



(£o+Po)tk 

The speed of diffusion of a mode with wavenumber k is now: 

d\w\ r 



v diffcausal (k) 

which does not violate causality: 



V 



dk 



(47) 



(48) 



lim Vdiff(k) 



V 



k— >oc 



(49) 



/ (£0 +Po)Tn 

where r n ^ 0. Beyond the Maxwell-Cattaneo law, the more complete formulation of vis- 
cous hydrodynamics is the Miiller-Israel-Stewart theory [12J, which contains the Maxwell- 
Cattaneo law as a limit. A more precise study is found in [T2] . Here we just wish to point 
out that there are alternative approaches, in which causality is preserved. 

We present a different approximation scheme, which goes beyond the linear approxima- 
tion, preserves nonlinear terms and does not violate causality: the reductive perturbation 
method, a technique that will be presented in the next section. When this method is ap- 



plied to the relativistic Navier-Stokes equation (33) we arrive at the conclusion that, for the 



purpose of studying perturbations, the equation: 



dv 
di 



+ [v- V)u 



(e + p) 



Vp + v 



dp 
di 



+ 



[e + p) 



,,V : -'r+ {( + - r] \v(V -if) 



(50) 



is equivalent to (33). In other words, considering (50) or (33) leads to the same wave 
equation. 
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E. Continuity equations 



1. Entropy 



The continuity equation for the entropy density is given by (25): 



(51) 



Using (§, Q, = and = 1 we can rewrite the above equation in the form: 



(52) 



Using (16) and u M = (7,7??) the last equation can also be rewritten as 

2 



ds ^± $7 ^± ^± _ r> (d^i 

7^ + 7Vs • v + s— + s V7 • v + 7s V • v = -— I — 



>1 



7 ' Wv lv ' 



|(9V)^ + ^(^ + c 



c?7 -> _ -* J 2 
— + 7V • v + V7 • v 
at 



(53) 



which is the relativistic version of the continuity equation for the entropy density s. In the 
case of an ideal fluid (rj = ( = 0) we recover the entropy density conservation: 



ds -* _ dj -* -* _ 

7— + 7 v s • v + s — — h s V 7 • v + 7sV • v = 
at at 



which, with the use of: 



becomes [10] : 



d"f 



dv 



dt ^ V m 



and V7 = 7 3 z;Vf 



+ 7 2 ^s i — + v - Vvj + V ■ (sv) = 



(54) 



(55) 



5. Baryon density 



The relativistic version of the continuity equation for the baryon density ps in ideal 
relativistic hydrodynamics is: 

d„jB V = (56) 



Since jb v = u h ' p B and using (54) the last equation can be rewritten as [9]: 



B 



dp 
~dt 



( 'dv - \ 
+ l 2 vp B [-^+v-Vv\ +V -(pbv) = 



(57) 
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F. Non-relativistic limit 



In what follows we recover the non-relativistic limit of the continuity equation for the 
entropy density, continuity equation for the baryon density and also for the Navier-Stokes 
equation. The non-relativistic limit is essentially given by v 2 << 1 and 7=1. 

1. Continuity equation 



Using v « 1 and 7 = 1 in the equation for entropy density (53) we find: 



f f + V • (sv) = -^(dV)d jUi + i (Jr? + () (V ■ v) 2 (58) 

For a perfect fluid r\ = ( = we obtain the usual continuity equation [9]: 

ds -> 

For the baryon density continuity equation we have, after applying v 2 « 1 and 7 = 1 in 
©: 

^ + V • (p B v) = (59) 

2. Navier-Stokes equation 

The energy density of the fluid element of mass M and momentum P is: 

E VM 2 + P 2 
e ~Vol~ Vol 

In the non-relativistic limit: M » P we have: 

M 



Vol Vol 

The volumetric density of fluid matter is given by ^ = p and hence: 

£ S p (60) 

The pressure is defined as the ratio of the force per area of the fluid element A. The force 
is given by the time derivative of the momentum P. The pressure is then: 

force 1 dP . 

v = ^ = -ATt (61) 
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and so 



For v 2 << 1 we find: 



and consequently: 



- ^dp l^dP 1 _d 2 P 

vv + v — = — V 1 — v 

p dt A dt A dt 2 



,->pp 










dt 



- _<9p 1 -idP 



With the use of (61) we finally obtain: 



-> dv -> 
Vp + v-£- = Vp 
dt 



(62) 



According to these last approximations and using (60) we conclude that e+p takes the form: 



e + p ^ p + 



1 dP 
A~dt 



and with P » ^^jr (because M » P) we find: 



e + p ^ p 



(63) 



Inserting 7 = 1, (62) and (63) in (33) we find: 

'dv 



P 



dt 



+ (v-V)v 



+ Vp 



V 2 v- V{V-v) 



dv ._ ^,.dv 



- ^ dv _ -* 
-V-u — -v- V 
at 



+ (C-g»7)^ -V(V-tT) 



• V)v 

dv 



s v '"-" v '« 



({7- V)u 



V-^ =0 



(64) 



Neglecting the terms of order v 3 , which are v ■ V (u • V)v 

'dv 



P 



dt 



and 
1 \ 



• V)v 



V-tJwe have: 



-(C + yi)|[^-s) 



Vp - rjV 2 v - [C+ V( V • 

d r, 



7/ 



at 



(65) 



We still need to estimate the relative size of the terms in last equation. This is made by the 



comparison of terms without viscosity to others with viscosity. From (65): 

77 d - 



- w 9 r _ - 
p(v-V)v — T] — (v-V)v 



M 
' Vol 
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y • v)?r — 



M 2 at 



(P- V)P 



where we have used the volumetric density of the fluid matter A^ = p and the non-relativistic 
momentum P = Mv for the element with mass M for the fluid matter. In the non-relativistic 
limit we have, as usual: M » P and the last expression becomes: 



-> d - — 

p{v-V)v — rj— (v-V)v = p{v-V)v 



(66) 



since 



M V d 

—(v ■ V)v » 



M 2 dt 



(P-V)P 



Analogously we also have from (65): 



dv 
7 dt 



d r 



since 



C+ 3 V Wt [v(V.v) 



1 d 



= P 



dv 
dt 



(67) 



M dv _( 1 
Voldt >> ~V + 3 V J M 2 8t 



P (V ■ P) 



The non-relativistic Navier-Stokes equation is then given by using (66) and (67) in (65) 



dv 

at 



+ (v- V)v 



X Vp+ -V 2 v+ 1 ((+ ^)V(V -v) 
p p p\ 3 ' 



(68) 



Again, the perfect fluid is described by r] = ( = in (68) and gives the non-relativistic 
version of Euler equation [10] : 

dv -> I -> 

— + (v-V)v = — Vp (69) 
ot p 



III. THE REDUCTIVE PERTURBATION METHOD 



A. Linearization 



We start from the equations of ideal relativistic hydrodynamics and, using the lineariza- 
tion approximation, we derive a wave equation for perturbations in the pressure. This 
equation has traveling wave solutions which represent acoustic waves. In the derivation pre- 
sented here we follow closely the reference [TT] . In the presence of perturbations the energy 



density and pressure for the relativistic fluid are written as (35): 



and 



e(f,t) = e + 5e(f,t) 



p(r, t) = po + 6p(f, t) 



(70) 



(71) 
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respectively. The uniform relativistic fluid is defined by Eq and po, while Se and Sp correspond 
to perturbations in this fluid. Energy-momentum conservation implies that for an ideal fluid 
IP 1 ' = and from ^ we have: 

dfjT^ = (72) 
where T^ u is the energy-momentum tensor given by: 



= (e + p)u»u v - pg 



(73) 



Linearization consists in keeping only first order terms such as Se, SP and v and neglect 
terms proportional to: 

v 2 , vSe, vSP, v-Vv, (v-V)v (74) 
and also neglect higher powers of these products or other combinations of them. Naturally 



we have 7 ~ 1 . From ( 72 ) we have 



u»d u [(e + p)u v ] + {e + p)u v d v u» - d„(pg Ufl ) = 
The temporal component (p, = 0) of the above equation is given by: 

7<9 [(e + p)i] + idi[{e + p)u 1 } + (e + p)u a d i + (e + p)u l da - d p = 



(75) 



(76) 



which, after using (74) and 7 ~ 1, becomes: 



d o {e + p) + d i [{e + p)v i ]-d o p = 



or 



de 
di 



V • [{e + p)v] = 



For the j-th spatial component (p = j) in (75) we have 



(77) 



u 3 d [{e + p)u u ] + u ] di[{e + p)u 1 } + (e + p)u d u j + (e + p)u l diV? - &>p = 



which, with the use of (74) and 7 ~ 1, becomes: 





Of 



[(s + p)v\ + Vp = 



(78) 



Substituting the expansions (70) and (71) in (77) and (78) we find: 





[e + Se] + V ■ [{e + de + p + 8p)v\ = 



(79) 
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and 



0_ 
dt 



[(e + Se + po + Sp)v\ + V[ Po + 5p] = 



Neglecting the terms listed in (74) in (|79|) and (80) these equations become: 



d(8e) 
dt 



and 



dv -> 
(e +Po)^ + V(Sp)=0 



(80) 



(81) 



(82) 



Equation (81) expresses energy conservation and equation (82) is Newton's second law. 



Integrating (82) with respect to the time and setting the integration constant to zero we 
find: 



V(Sp)dt 



which inserted in (81) yields: 



d(Se) 
dt 



(e + Po) 

V 2 (Sp)dt = 



Performing the time derivative we obtain: 

d 2 (Se) 
dt 2 



- v 2 (5p) = 



Assuming that 



de 

Se = —Sp 
op 



with de/dp being a constant, we have ([85]) rewritten as: 



- V 2 {Sp) = 



de d 2 (Sp) 



(83) 



(84) 



(85) 



(86) 



(87) 



dp dt 2 

The above expression is a wave equation from where we can identify the velocity of propa- 
gation as: 



de, 



where c s is the speed of sound. Equation (87) can then be finally written as: 

1 d 2 {Sp) 



V 2 (Sp) - 







(88) 



(89) 



c 2 dt 2 

which describes the propagation of a pressure wave in the fluid. This study ensures that the 
analysis of perturbations with the linearized relativistic hydrodynamics leads to the standard 
second order linear wave equations and their traveling wave solutions, such as acoustic waves 
in the hadronic medium. 
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B. Beyond linearization 



While linearization is justified in many cases, in others, where perturbations are not so 
small, it should be replaced by another technique to treat perturbations keeping the nonlin- 
earities of the theory. This is where a physical theory, in our case relativistic hydrodynamics, 
may benefit from developments in applied mathematics. Indeed, since long ago there is a 
technique which preserves nonlinearities in the derivation of the differential equations which 
govern the evolution of perturbations. This is the reductive perturbation method (RPM) 



We start the RPM description by considering the simple linear wave equation: 

dF a dF _ Q 
dt dx 



(90) 



which describes one-dimensional waves in an ideal fluid. In the above equation a is a 
constant. Let us expand of F(x,t) around the constant value F in terms of the small 
expansion parameter a (0 < a < 1) : 



F(x,t) = F Q + <jF x (x,t) + a 2 F 2 (x,t) + a 3 F 3 (x,t) + 



(91) 



Inserting (91) into (90) we obtain a series: 



J dF, Q 
dt dx I 



-^r n 



for j = 1,2,3,... 



(92) 



The coefficients of each power of a must vanish independently, i.e., each bracket must vanish. 
This condition yields a set of differential equations for the Ff 



3 + a^r 1 = 



(93) 



dt dx 

Now we consider the simplest nonlinear wave equation, the so called breaking wave equation: 



dF ^ p dF _ 
dt dx 



(94) 



where a is a real coefficient. Performing the expansion (91) the above equation becomes: 

(95) 



dF x 2 dF 2 aF\ 2 p dF 2 2 dF 1 

a—- + a — — + craF — h a aF — h a aF 1 — h 

at dt dx dx dx 







Let us look for simple nonlinear differential equations for the perturbations F±, F 2 , ... . We 
expect to find algebraic structures similar to (94) in equation (95) in order a and a 2 . To 
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do so, we introduce the new coordinates r and £, which are connected to t and x. This 
coordinate transformation is such that: 











a 



and 











a 



(96) 



dt dr dx <9£ 

where n and m will be determined. We note that a highly nontrivial aspect of this transfor- 



mation is that it contains the same small parameter a used in the expansion (]95|). Inserting 

dFi 



(96) into (95) we find the following equation: 



or 



+ a m+2 aF 1 -± + aaF 
or a£ 



m+l 







(97) 



where, for the sake of simplicity, we have chosen F = and we have truncated the sum 
keeping only the lowest order terms in a. We observe that if we choose: 



n + 1 = m + 2 



and 



n + 2 > 3 



(98) 



the first and third terms can be grouped together and, for the particular choice n = 3/2 and 

dFi 



m = 1/2, equation ([97]) becomes: 



a' 



<9F 2 
<9r 







(99) 



where, just before arriving at the above expression, we have divided the whole equation by 
a 1 / 2 . From the term proportional to a 2 we have: 



dFi „ dFi 

+ aFi^- = 



(100) 



Or ' «9£ 

as expected. Including all the higher order terms, F 2 , F 3 , ... etc, we would obtain a set of 



differential equations. Solving them we would be able to write the complete series (91). In 



practice, since a is small, this series is dominated by the first terms and it is often sufficient 
to compute only F\. 



From this simple exercise we conclude that for a nonlinear wave equation such as (94) it 



is necessary to introduce the "stretching" operators (96): 

d 



a 



3/2 







and also 



dt dr 

which are related to the "stretched" coordinates: 



d_ 
dx 



a 



1/2 



d_ 



(101) 



i = o^x 



and 



r = a 6 'H 



(102) 
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If there was a dissipative term in (94), the wave equation would have the following form: 

d 2 F 



OF F dF 

dt dx 



v 



dx 2 



(103) 



which is called Burgers equation with a dissipative (viscous) coefficient v. It becomes, after 
performing the expansion (91) and keeping only terms up to order a 2 : 

d 2 F l 



a 



dF 1 



a 



,dF2 
dt 



aaFr 



8F 1 2 „ dF 2 „ 
ox ox 



dFt 
dx 



vo- 



dx 2 



+ pa 



,d 2 F 2 
dx 2 



(104) 



As before, we try to find in equation (104) a structure similar to ( 103). Keeping the previous 



choice n = 3/2 , m = 1/2 and applying (101) to (104) we find: 



a 



5/2 



dFt „ dFt 
+ aFt 



dr 



0$ 



7/2 



dF 2 
dr 



2^1 , ...3 



+ va 



d 2 F 2 



(105) 



When dissipative effects are included we perform the following transformation in the dissi- 
pation coefficient [44"1 135] : 

u = a l/2 v (106) 



in (|105|), resulting in: 

dFt 



a 



5/2 



I 771 ^Ft 

d7 + aFl W 



+ a 



7/2 



dF 2 
dr 



vo 



d 2 F 2 



As before, we divide the last equation by cr 1 / 2 , obtaining: 

~dF 2 



a' 



dFt „ dFt 
+ aFt 



dr 



dt 



+ CT U 



dr 



~ 2^ Ft „ 3 

Per 2 —— + VO 6 



d 2 F 2 



The terms proportional to a 3 in (108) are neglected and from the a 2 terms we have: 

dFt , dFt ~d 2 Ft 



dr 



(107) 



(108) 



(109) 



which is a Burgers equation like (103). In addition, we conclude that when a dissipative 



term is present in a nonlinear wave equation, apart from (101) and (102), we also need the 



transformation (106) for the dissipation coefficient. 

If a more complete wave equation is considered, i.e., when nonlinear, dissipative and 
dispersive terms are present, we have the Korteweg-de Vries Burgers (KdV-B) equation: 

d 2 F 



dF +a F dF +(3 d3F 
dt dx dx 3 



v 



dx 2 



(110) 



Now, (3 is the dispersive coefficient. Performing the same calculation we naturally find: 



a 



dFt „ dFt n d 3 Ft 



dr 



dt " de 



a' 



d 2 Ft 



111) 



20 



which gives the KdV-B: 

~d7 + aFl W + ^W = u W (112) 

When more dimensions and different coordinate systems are considered, the procedure can 
be systematically improved. 



C. Some special cases 



In the framework of the RPM we transport the equations of hydrodynamics from the space 
of cartesian, spherical or cylindrical coordinates to the space of the "stretched coordinates" . 
Some well known equations and the coordinate transformations and expansions required to 
obtain them are given below. 



1. One dimensional KdV equation 

We write this equation in cartesian coordinates (x), in radial cylindrical coordinate and 
in radial spherical coordinate (r). The corresponding "stretched coordinates" are given by 
£ for space and r for time as: 



and 



a 1 / 2 

Z=—(X-c s t) (113) 



^3/2 

r = — c s t (114) 



where X = x or X = r. The above transformation of coordinates must be made simultane- 
ously with the following expansions of the baryon density, energy density and fluid velocity 
around their equilibrium values: 

p= — = l + a Pl + a 2 p 2 + ... (115) 
Po 

e = — = I + ae 1 + a 2 e 2 + ... (116) 

v = — = CTVi + <T 2 t> 2 + . . . (117) 
C s 

For two and three dimensional nonlinear wave equations we will use the "stretched coordi- 
nates" presented in 43J and references therein. 
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2. Two dimensional cylindrical breaking wave equation 

In this case the coordinate transformation is given by: 

a 1 / 2 a a 3 ' 2 

R=—(r-c s t), Z=-z, T=—c s t (118) 

The energy density and fluid velocity components are expanded around their equilibrium 
values: 

e = — = l + ae 1 + a 2 € 2 + a 3 e 3 + ... (119) 

v r = — = av ri + a 2 v r2 + a 3 v r: + . . . (120) 
c s 

v z = ^ = a 3 ' 2 v zi + o h l 2 v Z2 + a 7 l 2 v z , + ... (121) 

C s 

3. Three dimensional cylindrical KP equation 
In this case we have: 

1/2 3/2 

R= a j-(r-c s t) , $ = o- x >\ , Z = ^z, T=°—c s t (122) 

and for the baryonic density and fluid velocity components expansions are: 

p = — = 1 + api + o~ 2 p 2 + o- 3 p 3 + ... (123) 
Po 

v r = — = av ri + a 2 v r2 + <j 3 v r3 + . . . (124) 

Cs 

= — = o 3 ' 2 v^ + a 5/ \ 2 + o- 7 /\ s + ... (125) 

Cs 

v z = V - = <j 3 ' 2 v zi + a" 2 v Z2 + a 7 ' 2 v z , + ... (126) 

Cs 

p 4 /3 = [1 + ( api + a *p 2 + . . ,)] 4 /3 ^ l + 4 (jpi + 4^ + (127) 
p = [1 + (<Tpi + (T 2 P2 + . . .)] 1/3 = 1 + ^Pl + ^ 2 P2 + • • • (128) 
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4- Three dimensional cartesian KP equation 



Here the stretched coordinates are: 



a 1 ' 2 , s a a a 3 ' 2 , , 

X = —(x-c a t), Y = ~y, Z = l z > T = ^ Cst ( 129 ) 



and the expansions are: 



p = — = 1 + °pi + a 2 p 2 + a 3 p 3 + ... (130) 
Po 

v x = — = crv Xl + o- 2 v X2 + a 3 v X3 + ... (131) 

C s 

Vy = — = (T V2 V yi + (J 2 V y2 + <X 5/ \ 3 + . . . (132) 
Vz r 3/2„, ,2, ,5/2 



C 



o ill v Zl + cr 2 v Z2 + a b/I v zz + ... (133) 



p 4/3 = l + ^ Pl + ^ 2 P 2 + ... (134) 
p 1/3 = l+3api + Vp 2 + ... (135) 

In all cases L is a characteristic length scale of the problem, e is the equilibrium (or 
reference) energy density, po is the equilibrium (or reference) baryon density and c s is the 



speed of sound. When viscosity is included, we perform the transformation (106) as in 
HUSH]: 

C = o- 1/2 C and also 77 = cx 1/2 f\ (136) 

Once the equations are in the "stretched spaces" and expanded, we neglect terms propor- 
tional to cr n for n > 2 and organize the equations as series in powers of er, a 3 / 2 and a 2 . 
These equations form a system of differential equations which are combined to yield the 
final nonlinear equation for the relevant perturbation. Finally we transform the nonlinear 
equation back to the cartesian (cylindrical or spherical) space and solve it. 



IV. THE EQUATION OF STATE 

The equations of hydrodynamics discussed in the previous sections must be supplemented 
with an equation of state (EOS), i.e., a relation between pressure (p) and energy density (e) 
or matter density (p). In relativistic hydrodynamics, we can write the EOS as: 

P = c 2 e (137) 
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where the speed of sound c s must be smaller than one. The equation of state is derived 
from microscopic theories of the strongly interacting system. As mentioned above, the 
fundamental theory of strong interactions (QCD) predicts that cold and/or dilute systems 
are in the hadronic phase, where the degrees of freedom are baryons (proton, neutron, A, ...) 
and mesons (it, p, ...). At higher densities and/or temperatures, there is a phase transition 
and the formation of the quark gluon plasma (QGP), a phase where these particles are 
free and all hadrons are dissolved. In all phases we may have thermal equilibrium and 
hydrodynamical behavior. In both hadron and quark-gluon fluids we may have excitations 
and the nonlinear propagation of perturbations. In what follows we shall study the formation 
of these nonlinear waves and the role played by the equation of state. 

In the pioneering study of Refs. [201 SO] the authors studied the equations of non- 
relativistic hydrodynamics and, with the help of the RPM, they were able to derive a KdV 
(Korteweg de-Vries) equation for a perturbation in the nuclear density. A key ingredient 
in that work was the equation of state, which established the following relation between 
pressure (p) and density (p): 

y = V0 (138) 
where is a potential, playing the role of heat function, given by: 

0= — L' + c 2 VV + -l (139) 
Po 1 J 

where c\ and C2 are constants and p' is the deviation of the density from its equilibrium value 



p' = p — po- The Laplacian in (139), combined with the gradient in (138) gives origin to 
the cubic derivative of the KdV equation. As a result the authors arrived at the conclusion 
that a KdV soliton may exist in cold nuclear matter, when, for example, a light nucleus 



is impinged on a heavy nucleus. This conclusion relied strongly on (139) and (138), which 
come from an oversimplified description of the nuclear interactions. In [35H38] the nuclear 
interactions were described in terms of a relativistic mean field model, which is a variant of 
the Walecka model. Within this framework we may try to answer the question: what is the 
microscopic origin of higher order derivative terms appearing in the equation of state? This 
will be discussed in the following sections. 
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A. Hadronic Matter 



In this subsection we present an equation of state derived from a successful relativistic 
mean field model: Quantum Hadrodynamics (QHD) or Walecka Model. For a modern 
approach, using chiral power counting in an effective field theory for nuclear matter with 
nucleons and pions as degrees of freedom, see Ref. [4?] . 

This model is well established and is the subject of textbooks as Refs. |3j . The Lagrangian 
density of nonlinear QHD is given by: 

C = jky^id* - g v Vn — (M — g s ^ + ~(d^cf> - m s 2 (f> 2 ) - ~ V"+ 

+ \m v 2 V^ - b -<p 3 - + C d (140) 
where = <9 M K - d v V^ and: 

C d = d^-Md^V^ (141) 
my 

Except for the last term, this is the standard nonlinear QHD Lagrangian which is able to 
reproduce all the main features of nuclear matter and finite nuclei. This last term was added 
only in [35H38] and illustrates how to include higher order derivative terms in the equation 
of state. 



The Lagrangian (141) is the modern version of (139) and (138). In (140) the degrees 



of freedom are the baryon field ip, the neutral scalar meson field and the neutral vector 



meson field V^, with the respective couplings and masses. The last and new term (141) is 
designed to be small in comparison with the main baryon-vector meson interaction term 
gviplfiV^ip- Because of the derivatives, it is of the order of: 



* ~ 0.12 (142) 

mi/ mi 



by liuy 

where the Fermi momentum is hp ~ 0.28 GeV and my ~ 0.8 GeV. The form chosen for 
the new interaction term is not dictated by any symmetry argument, has no other deep 
justification and is just one possible interaction term among many others. It is used here 
as a prototype to study the effects of higher derivative terms, which, as it will be seen, may 
generate more complex wave equations, such as KdV. The parameter d is free and plays the 
role of a "marker". Setting d equal to zero switches off the new term and we recover the 
usual QHD. On the other hand d — 1 means that the coupling gy is the standard one. Other 
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values imply a correction in this coupling. As mentioned in the beginning of this section, 
interesting phenomena, such as KdV solitons, may appear as a consequence of the use of 
equations of state with higher order derivative terms. 

Baryon number propagation in nuclear matter satisfies the diffusion equation: 

dps 



dt 



(143) 



where the diffusion constant D has been numerically calculated and studied as a function of 
density and temperature. For example, in jl8] it was found that D ~ 0.35 fm at densities 
comparable to the equilibrium nuclear density and temperatures of the order of 80 MeV. 
This number is small compared to any nuclear size scale and can be interpreted as indicating 
that 

rln r, 

(144) 



« V p B 



dt 

and therefore the density gradients do not disappear very rapidly in nuclear matter. Because 



of the above inequality we can neglect the time derivatives in (141) 



The usual mean field theory (MFT) approximation is based on two assumptions: a) the 
baryonic sources are intense and their coupling to the meson fields is strong and b) infinite 
nuclear matter is static, homogeneous in space and isotropic. 

The first assumption above is implemented performing the substitutions: 



and 



<f) -K <f) >= 0o 



(145) 



(146) 



in 



(140) and obtaining C* = C{y^ -K >= S^Vo <f) ~*< 4> >= 4>o)'- 

C* = i>[(i^d» - g vlo V ) - (M - gsfoM + Ud^ d^ - m s 2 4> 2 ) + 



+ ^Vo) 2 + ^m v 2 V Q 2 - -0o 3 - |0o 4 



The equations of motion are given by [38 , |4*9H5T] : 

dC* dC* 

— _ Q —— + d u d 

where r\i = (f) Q , Vq and read: 



dC* 



d{d„d v r)i) 



- VV Q + m v 2 V = g v ipl ^ 



(147) 



(148) 



(149) 
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{d^ + m s 2 )4> = ggifa) - 60 o 2 - 



i^nd* - gvloVo — (M — g S (fio) 



* = 



(150) 
(151) 



The effective mass of the nucleoli is given by M* = M— g$ 4>o- In the case of pure QHD ( 147) 



because of the interaction between the nucleon and the vector meson, we can anticipate that 
the second order derivative term in the field V^, which comes from its equation of motion, 
may be transferred to the term ip'-f^ip and hence to the baryon density ps- In the mean field 
approximation, only the time component of the field, V°, contributes, yielding the term 
V 2 Vq. It is possible to estimate V 2 Vq from the equation of motion (149). To do so, we first 



rewrite the equation ( 149 ) as function of the baryon density: 



VV + m v 2 V = g v p B 



(152) 



Next we assume that V 2 Vo/m 2 / << Vo, neglecting the Laplacian in (152) to find the first 
order estimate of V^: 

(153) 



rr 9V 

Vo = -^Pb 



my 1 



An improvement of this estimate of Vq is obtained taking the Laplacian of (153) and substi 



tuting it in the Laplacian present in (152). After this simple algebraic procedure, we solve 



the resulting equation for V and obtain: 

Vo 



9v 9v ^2 

ZpB + 1 V pB 



my Vtly 

The energy- momentum tensor is given by [38l I49TI5T] : 

dC* 



dC* 

T^ = —-—(d v Vi )-g^C* 
The energy density is: 



which becomes 138 



0, 



didf.dprji) 



dC* 



e =< T 00 > 



(154) 



(155) 



(156) 



e = 2 Wo) 2 + I (V0 O ) 2 + \gv PB Vo + \rnsV + ^ + \^ + j^y 3 fj <Pk(J? + M* 2 ) 1 ^ 2 

(157) 



Inserting (154) in the expression of the energy density (157) and also writing 0o m terms of 



M* we find the following expression: 



9v 2 , 9v z 



2m v 2 ' 



2m y 4 ' 



— W + 7^— d PBV 2 p B + ^(M- M*) 2 

Is 



2gs 2 
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(M-M*) 3 (M — M*) 4 
& ~ n r c - -. r 



d 3 k(k 2 + M* 2 ) 1 ' 2 



(158) 



iS- W ' (2vr) 3 Jo 

In a first approximation, the variables ps, V 2 p# and M* are independent from each other 
and therefore, taking the derivative of the above expression with respect to we have: 



de g v 2 gv 2 ^ 2 



(159) 



dps my 1 - £Uly 

We will add the two sources of inhomogeneities in ps (which are responsible for a non- 
vanishing V 2 ps)- For cold nuclear matter we have then [38J: 



If d 
2\di 



(M-M*)!] 2 I 

gs 



(M-M*)]) 2 , m s 2 



gs 



+ ^-(M-M*) 2 + 



(M-M*) 3 (M — M*) 4 
b- — - — ir J - + c- 



^gs A 



+ 



is 

Is 



gv' 

2m 



;Pb 2 +[d + 



v 



2 J m v ' 



PsV p B + 



kp 



d 3 k(k 2 + M* 2 ) 1 / 2 



(160) 



(2tt) 3 Jo 

where 7 S = 4 is the nucleon degeneracy factor. For hot nuclear matter, the energy density 
as described in [38J is given by: 

1 f d 



2\dt 

(M-M 
'~ 3^ 



(M-M*)^ 2 1 
+ 2 I 



gs 



(M-M*)]) 2 , m s 2 



gs 



+ — > (M-M*) 2 + 

2gs 2 



*\3 



+ C 



+ 



(M-M 



*\4 



+ 



gv' 

2my 2 



Pb +[d + 



2) my' 



PbV 2 pb+ 



Is 



(27T) 



J d 3 k h + [n n {T, v) + n % {T, v)\ 



The baryon density is: 



Pb 



Is 



with 



(2vr) 3 

ng(T, v) = 
nAT, v) = 



d 3 k [ ni (T,u)-n n (T,u)] 



1 



1 _|_ e {h+-u)/T 
1 

1 _|_ e {h + +v)/T 

gv 



and 



u = hb- g v v + d-^(d*d v v ) 

my z 



h + = (k 2 + M* 2 f/ 2 



(161) 

(162) 

(163) 
(164) 
(165) 

(166) 
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It is important to note that the term Cd provides the "(d+l/2)-term" in the energy densities 



(161) and (160) that generates the KdV equation. The nucleon effective mass (M*) at zero 
temperature is obtained through the minimization of e with respect to M*: 

d£ 



DM* 



which, with the help of (160) yields: 



M* = M - 



9s 7s 
m s 2 (27r) 3 Jo 



d 3 k- 



M* 



(fc 2 + M* 2 )V2 



+ 



9s' 
ms 2 



(M - M*f + —AM - M*f 



9s A 



9s H 



Analogously, at finite temperature: 



M* = M 



9s Is 
m s 2 (2vr) 3 



d 3 k 



M* 



n^T, v) + n n {T, v) 



+ 



9s 



2 r b 



*\3 



(M - M*y + — -(M - M*) 



(167) 



(168) 



-9s° 9s H 

and we conclude by these two last expressions that "cf-term" does not affect the nucleon 
effective mass. 



B. Quark Gluon Plasma 

The idea that quarks and gluons may exist as free particles in a deconfined phase was 
advanced long ago [H2] in the context of compact stars. In these objects gravitation is strong 
enough to compress the nucleons and make them overlap with each other. The distance 
between the quarks can be so small that, due to asymptotic freedom, they almost do not 
interact. In this picture, we may have large regions of space populated with free quarks. 
This state is called the (cold) quark gluon plasma. This naive description of QCD in extreme 
conditions of density was later revisited with much more sophisticated models. Quarks and 
gluons may also form a hot quark gluon plasma. The hot QGP has been extensively studied 
in lattice simulations and also in heavy ion experiments at CERN and LHC Today we know 
that a hot and deconfined phase is formed in the existing accelerators, but it is far more 
complicated than previously imagined. In particular the quarks and gluons are not really 
free. Instead, they still interact strongly with each other forming a state called the strongly 
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interacting QGP, or sQGP. Also the non-trivial vacuum structure persists until relatively 
large temperatures. 

A starting point to study in a unified way both the cold and hot QGP is the MIT bag 
model [53]. According to this model in its simplest version, massless and non-interacting 
quarks live in a spherical cavity (the "bag" ) in the physical QCD vacuum, which is a medium. 
The confining property is represented by a constant term, "B", called the bag constant. In 
the next subsection we shall briefly show how to calculate the QGP equation of state with 
the help of this model. 

1. The MIT Bag Model 

In what follows we consider only quarks u and d. Each quark has three color states and 
they are massless. We also have eight massless gluons and we neglect the interactions in 
the QGP. The confinement property is included in the model through the introduction of a 
constant, positive energy per unit volume in the vacuum: 



that can be interpreted as the energy needed to create a bubble or bag in the vacuum, in 
which the noninteracting quarks and gluons are confined. B is known as "bag constant". 

The parameter B can be extracted from a phenomenological analysis of hadron spec- 
troscopy or from lattice QCD calculations. There is a relationship between B and the 
critical temperature of quark-hadron transition T c which is determined by considering that 
during the phase transition the pressure is zero. 

The quarks have baryon number 1/3 and the chemical potential for gluons is zero. The 
baryon density, energy density and pressure for the QGP are given by: 




(169) 




(170) 



where 



n k = n k( T ) 



1 



(171) 



1 _|_ e (k-±fj, B )/T 



and 



H = %( T ) 



l 



(172) 



1 _|_ e (k+^ B )/T 
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where hb is the baryon chemical potential. The energy density is given by: 

e = B + ^ / d 3 k k (e k / T - l)- 1 + j^Jd 3 kk [n s + %] (173) 

where the first term of the expression above is the gluon contribution. The pressure is given 
by: 

(174) 



P 



3 {(2ir^ 
The degeneracy factors are: 



7g = 2 (polarizations) x 8(colors) = 16 for gluons 



(175) 



and 

7q = 2(spins) x 2(flavours) x 3(colors) = 12 for quarks (176) 

The integral of the gluon distribution function can be calculated analytically and the ther- 
modynamics of QGP can be summarized in the following expressions: 

2 roo 



and 



I'l> = ^2 J Q k [rag - n- h 



87r 2 6 r°° 

3(p + B) = e-B= — T 4 + — / dk k 3 [nz + n 
15 n 2 Jo ft 



(177) 



(178) 



]_5 ' -2 /,. ' ' 1 feJ 

which provide us with the EOS of QGP with all baryon densities and at all temperatures: 



1 4 « 

p = -e — -B 



(179) 



The sound speed c s is given by: 



dp 

de 



180) 



2. The cold QGP 



In the core of a neutron star, the temperature is zero and baryon density is considerable. 
This is a good place to study baryon density perturbations. The quark distribution function 



becomes the step function with \xb = 3kp and (177) becomes: 



Is 

6tt 2 



Pb = —k F 3 



(181) 
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Expression (178) becomes: 



3(p + B) = e - B 



6 (k 



7T 2 V4 



k p 



2tt' 



:k F 4 



;i82) 



Using (181) we find 



3(p + B)= £ -8=Q 7/ VV 3 



and so 



e(PB)=Q 7/ \ 2/ W /3 + B 



;i83) 



Inserting (183) into (179) we find: 



p(PB) = l(l) V3 n 2/ W /3 -B 



184) 



From ( 183 ) and ( 184 ) is possible to rewrite the sum of energy density and pressure as: 



From (179) we have: 



Vp = -Ve 



and also 



dp 1 de 
di = 3di 



Using (183) in (186) we find: 



Vp 



B 



and 



dP = V3\^ 2/3 1/3 dps 
dt 9\2J PB 8t 



;i85) 



;i86) 



(187) 
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3. The hot QGP 



The hot QGP is formed in heavy ion collisions in the central rapidity region where we 
have pb = 0. For our purposes the most relevant physical quantity is the energy density 
e. Since p# = 0, the baryon chemical potential is zero ((jlb = 0) and so the distribution 



functions given by (171) and (172) are the same: 



n k- n k- ]_ _|_ e k/T 



and then (178) takes the form: 



87r 2 12 r°° k 3 

3(p + B)=e-B=^-T' + -J dk M fc . ._. 

15 7T 2 Jo [1 + e fc / T ] 



189) 
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We can perform the above integral analytically to find: 
From basic thermodynamics we know that 

-(£), 



(190) 



(191) 



which, with the use of (190) leads to the specific form for s: 



8 = ( ~ B + -^T 4 ) = 4 37 7r 2 T 3 
dT V 90 / 90 



(192) 



The parameter £>, the "bag constant" can be defined at the temperature Tg. In the bag 
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surface, (190) is given by: 



B = —ii 2 (T B ) A 
90 V ; 



Choosing £ 1/4 = 170MeV corresponds to T B = 91 MeV. Rewriting ( |190| ) as: 

.1/4 



T 



30 
37^ 



(s-B) 



(193) 



(194) 



and inserting it into (192) we find: 



/ \ ,37 2 



30 
37^ 



n3/4 



(e-B) 



In a compact notation: 



where 



s(e) = A(e - Bf /A 

3/4 



A = 4 — 7T 2 

90 



30 
37^ 



(195) 



(196) 



C. Mean field theory for Quantum Chromodynamics 

In spite of its phenomenological success the MIT bag model gives a poor representation of 
the quark gluon plasma. From heavy ion collisions there is convincing evidence that quarks 
and gluons interact strongly forming rather a "prefect fluid" than an ideal gas. Therefore, 
the picture of free partons must be modified. There are several ways to do that. Here we 
discuss the approach proposed in [41J, which can be called mean field QCD (MFQCD). It 
allows us to start from QCD Lagrangian and derive an equation of state, which incorporates 
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the effects of a sizeable strong coupling constant and also residual non -perturbative effects 
from the QCD vacuum. As an interesting by-product this equation of state supports the 
existence of KdV solitons in a cold QGP. In what follows we summarize the basic ideas 
of this approach and derive the expressions for the pressure and energy density, which are 
relevant for the hydrodynamical study. 

We first introduce the mean field approximation for QCD, extending previous works along 
the same line [54"| [55] . We consider a system of quarks and gluons which are represented by 
the QCD Lagrangian density: 

1 Nf r 

Cqcd = — a F; v F^ Ij;^ ^ - i<fl*G$ - S ijmq U] (197) 

4 9=1 

with 

F ap, = QHQav _ QVQa^ + g jabcQb^Qcu (198) 

where ipf and Cr° represent the quark and gluon fields respectively. The summation on q runs 
over all quark flavors, m q is the mass of the quark of flavor q, i and j are the color indices 
of the quarks, T a are the SU(3) generators and f abc are the SU(3) antisymmetric structure 
constants. For simplicity we consider only light quarks with the same mass m. Moreover, we 
drop the summation and consider only one flavor. At the end of our calculation the number 
of flavors will be recovered. Following [5H |55] , we shall start writing the gluon field as: 

Qan = A an + a a„ ^gg) 

where A afl and a aM are the low ("soft") and high ("hard") momentum components of the 
gluon field respectively. We will assume that A afl represents the soft modes which populate 
the vacuum and a atJ- represents the modes for which the running coupling constant is small. 
In a cold quark gluon plasma the density is much larger than the ordinary nuclear matter 
density. These high densities imply a very large number of sources of the gluon field. Assu- 
ming that the coupling constant is not very small, the existence of intense sources implies that 
the bosonic fields tend to have large occupation numbers at all energy levels, and therefore 
they can be treated as classical fields. This is the famous approximation for bosonic fields 
used in relativistic mean field models of nuclear matter [3J. It has been applied to QCD in 
the past and amounts to assume that the "hard" gluon field, a" is simply a function of the 
coordinates [3j: 

a a ^(x,t) =6^a^(x,t) (200) 
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with d u a a ^ ^ 0. This space and time dependence goes beyond the standard mean field 
approximation, where a a is constant in space and time and consequently d v a a = 0. We 
keep assuming, as in [H] , that the soft gluon field A atl is independent of position and time 
and thus d u A afi = . Following the same steps of |41j we obtain the following effective 
Lagrangian: 



C 



m G 



■-a a (V 2 a a )+ 2 



-OlnOl 



0"0 



B 



QCD 



(201) 



where the constant Bqcd is the same bag constant for QCD as defined in jUJ. In fact, 



the effective Lagrangian (201) is quite similar to the effective Lagrangian obtained in j4TJ. 



The new feature of (201) is the term -^^(V a§). For simplicity, we take the quarks to 



be massless. From the above Lagrangian, using (155), it is straightforward to derive the 



energy-momentum tensor T„ v , which gives us the energy density e and pressure p: 



■-^W + (S)»*» + (iS)**™»> 



( 27g 2 
\16m G 6 



V 2 p B V 2 (V 2 p B ) + Bqcd + 3 



1Q k F ' 
2tt 2 4 



(202) 



and the pressure is: 

( 27s 2 



V 



\16mG 2 



Pb 



9g 2 



( V 



\16m. 



G 



Pb V p B - 
V 2 p B V 2 p B 



9g 2 



Pb V 2 (V 2 p B ) 
V 2 p B V 2 (V 2 p B ) 



9g 2 



&QCD + — 



4 



m G 
9g 2 

Vp s • V(V 2 p B ) 



2. 2 4 C 203 ) 
where 7q is the quark degeneracy factor 7q = 2 (spin) x 3 (flavor) = 6 and hp is the Fermi 
momentum defined by the baryon number density by ps = kp' 3 /n 2 . The other parameters 
g, m G and Bqcd are the coupling of the hard gluons, the dynamical gluon mass and the 
bag constant in terms of the gluon condensate, respectively. An improved version of the 



EOS of (202) and (203) was used in the study of three dimensional solitons in cold QGP. 



These solitons are solutions of the Kadomtsev-Petviashvili (KP) equation, which is the three 
dimensional generalization of the KdV equation. The complete description of the calculation 
may be found in |42j . 
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V. NONLINEAR WAVE EQUATIONS 



In the previous sections we presented a review of the equations of hydrodynamics and we 
introduced the equations of state, which represent the microscopic dynamics of the corre- 
sponding fluids. Furthermore, we have presented a mathematical prescription (the RPM) to 
study perturbations in these fluids preserving the nonlinearities of the differential equations. 
The application of the RPM to several systems of interest leads to the nonlinear differential 
equations which we discuss in this section. Most of these wave equations were developed in 

[HSHHS M EES SSJ. 



A. KdV equation in nuclear matter 



Cold nuclear matter 



We insert (160) in the Euler equation given by (34). We then combine this Euler equation 
with (57) and follow the RPM procedure to obtain the following KdV equation: 

dpi. d Pi,fn 2n ~ d Pi 9v 2 po \d 3 pi pi , . 

where Q is a "geometrical factor" : 

, for cartesian coordinates: X = x 
Q e= { (205) 
1 for spherical coordinates: X = r 



The non-relativistic limit is obtained by the approximation 3 — c s =3 and so (204) becomes: 



dpi dp! - dp x ( l\( g v 2 p \d 3 p! pi 
which could have been obtained with the use of the same procedure applied to the nonrela- 



tivistic Euler equation (69) with p = MpB- 



2. Hot nuclear matter 



We follow the same steps listed above, but changing the (160) by the EOS (161) to find 
a KdV given by: 

dpi , dp\ ( 2 p B m v 2 c s 2 \ ^dp\ ( l\( c s \d 3 p\ p\ 
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In the last two wave equations p\ = api as denned in (115) and Q is denned by (205). 
Choosing d — — | in the above equations we eliminate the cubic derivative term and obtain 
a breaking wave equation. 



B. Breaking wave equation in QGP 



1. Cold QGP 



Inserting the cold QGP MIT EOS, ( 182 ) to ( 188 ), into the ideal hydrodynamical equations 



(|34|) and (|57|) and following the RPM procedure we obtain the following breaking wave 

dpi 



equation: 



dpi 
dt 



( _ + 2 c . dpi 
s dx 3 s dx 







(208) 



where again pi = api. 



2. Hot QGP 



Now we insert the hot QGP MIT EOS relations, (189) to (196), into the ideal hydrody- 



namical equations (34) and (55) and following the RPM procedure we obtain the following 



breaking wave equation: 



dii dii 
dt s dx 



1 + 



T 



In 



dii 



(209) 



where ii = aei. 



C. KP and KdV equations from MFQCD 



The energy density and pressure are given respectively by (202) and (203) and the wave 



equations are for perturbations in the baryon density given by pi = api as defined by (123) 



and (130). We use the ideal hydrodynamical equations (34) and (57) following the RPM 



procedure. In the calculations we find the following relation: 

/27^W 



(27gW 
V 8m G 2 



c s 2 + Svr 2 / 3 po 4/ V 



+ 7T 2 / V /3 



.4 



(210) 
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which relates the speed of sound to the background density po : 



c s 2 = -} ^ (211) 

+ 3^/3 po 4/3 



The wave equations are the following. 



1. Cylindrical KP equation 



dr 1 dt s dr 12 



7T 2 / 3 p 4 / 3 

3A 



$g 2 po 2 c s 



8m G A A 



d 3 Pi Pi 
dr 3 2t 



1 d 2 pi c s d 2 pi 



2c s t 2 dip 2 



(212) 



2. KP equation 



dx\ dt s dx 



a-cs 2 )- p0 



+ 



c s d 2 pi 
2 dy 2 



+ 



3A . 

c s d 2 Pi 
2 dz 2 



. d Pi . 



9ff 2 Po 2 c s 
8m G 4 A 



d 3 Pi 



(213) 



5. iMF equation 



The one dimensional cartesian particular case of (213) is obtained by neglecting the y 



and z dependence, so that (|213|) becomes the KdV : 
dpi 



+ c— + 
dt s dx 



(2-c 2 ) ( 27g 2 p 2 \ (2c s 2 -l) TT^Vg 
2 V 8m G 2 ) 2A 

9g 2 po 2 c s ] d 3 p 



A 



c 2 - 



+ 



8m G A A 



dx 3 







Taking the limit m G — > oo we obtain from (210) and (211): 

A = n 2 / 3 p ^ , 



1 

3 



and (214) becomes: 



— + c — -t 
dt s dx 3 



2 „ dpi 

CsPl 



dx 







c s pi 



dpi 
dx 



(214) 



(215) 



and we recover exactly the result (208), the breaking wave equation for pi at zero tempera- 



ture in the QGP derived from the MIT equation of state. 
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4- Breaking wave equation 



Neglecting the spatial derivatives in ( 202[ ) and (203), the equation (214) reduces to: 
(2-c s 2 ) (27g 2 p 2 \(2c s 2 - 1) tt'-W '/ , I 



dt s dx 



2A 



A v Cs 6 



c s p^ = (216) 



which is also a breaking wave equation for p\ with the po, and g dependence in its 
coefficients. 

D. KP-Burgers equation in hot QGP 



The relations (189) to (196) for the hot QGP are now inserted in the relativistic viscous 



hydrodynamical equations, i.e., the Navier-Stokes (33) and the continuity for the entropy 



density (53). The ideal case is recovered when the viscous coefficients are set to zero. 



Following the RPM procedure we obtain the following wave equations: 



1. Cylindrical KP-Burgers 



d_{dii 
dr \ dt 



dii 



or 



1 + 



Tb 
To, 



dei 
Or 



2t 2T \s 3s) dr 2 



1 d 2 ii c s d 2 ii 
+ 2cJ? Ikp^ 2 + ~2~d^ 



(217) 



2. Cylindrical Burgers 



Neglecting the ip and z dependence, the equation (217) becomes: 

i 



dix ggi c s 
dt +Cs dr + 2 



1 + 



T, 



B 



T, 







9il ' * 1 fC + t V -\ d ^ (218) 



£l dr + 2t 



2T \s 3 s dr 2 



Setting r] = ( = (ideal fluid) the (218) becomes: 



dii dii c s 
dt +Cs dr + 2 



1 + 



B 



To 



. dii . e\ n 
^ + 2t=° 



(219) 
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VI. ANALYTICAL SOLUTIONS OF NONLINEAR WAVE EQUATIONS 



We present some cases where particular analytical solutions exist. For the KdV equation 
we have the soliton solution. Soliton or solitary wave is a localized pulse which propagates 
without change in shape. For a detailed study of the KdV solitons we recommend the reading 
of [56]. The soliton will be also used as an initial condition in the study of the numerical 
solution of the spherical KdV and the breaking wave equations. For the Burgers, cKP, KP 
and cylindrical KP-Burgers we present the exact solutions as developed in [2U |25j [281 |30j [571 - 
[UO] where several techniques to solve nonlinear wave equations are presented. 



A. KdV equation in nuclear matter 



We first show the soliton solution of the KdV equations (204) and (206) in cold and hot 



nuclear matter, respectively [35T438] . We have only soliton solutions in the cartesian case 



Q = 0, so X = x. We also choose d = 1/2 as in [3E]- The equation (204) can be integrated 
and solved exactly and its soliton solution is given by: 

2 



pi{x,t) 



3(u — c, 



(3-c, 



2\-l 



sech 



my 
9v 



(u - c s )c s M ^ 
2p 



x — ut) 



and the solution of (206) is found by setting 3 — c s = 3 in (220): 

2 



pi(x,t) 



[u - c. 



-sechf 



9v 



(u - c s )c s M 
2p 



(x — ut) 



Finally, the solution of (207) is: 
3(u — c 



pi(x,t) 



2 n B m v 2 c 2 1 

L C s 



2gv 2 po 



sech 



[u — c s jmv 



2c, 



[x — ut) 



B. KP equations in cold QGP 



The cKP equation (212) has the exact analytical soliton solution [42] : 



pi(r,tp,z,t) = J^ sech2 \~^ 



ar + oz — \u + a - t 



(220) 



(221) 



(222) 



(223) 



where u is a parameter which satisfies u > ac s + b 2 c s /2a and the phase velocity given by 
nits appearing in 
u — ac s — b 2 c s /2a 



u + a^^- The constants appearing in the above expression are: 



hi 



a 3 c 2 



and 



Cl 

3a 2 c 2 



(224) 
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where 



and 



2 l s J 



c 2 = 



9jT Po c s 
8m G 4 A 



(225) 



(226) 



For the KP equation (213) we have the following soliton solution [32]: 



- i \ 3(?7 — w) )2 

pi[x,y,z,t) = — — — sec/i 



'(U-w) 
4A 3 c 2 



Ax + By + Cz- Ut 



(227) 



where A, B, C are real constants and w is given by: 



B 2 c< C 2 c< 



w = Ac s H — - — h 



(228) 



2 2 

We shall consider A > for simplicity and the parameter U such that U > w . For the KdV 



equation (214) we have: 



pi(x,t) = — — sech 2 

c 3 



[u - c. 



4c 4 



(x — ut) 



(229) 



where u is an arbitrary supersonic velocity and the constants C3 and C4 are given by: 



C3 



(2-c, 2 ) (27g 2 p Q 2 \{2c 2 -l) 7r 2 / 3 P o 4/ V 2 1 
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(230) 



and 



C4 



V Po c. 
rric^A 



(231) 



C. KP-Burgers equation in hot QGP 



The cKP-B (217) has the solutions: 



v 2DA ( C, 4r/ N 
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(232) 
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£i(r,z,(p,t) 
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(233) 



or the following solutions: 



x 2D A ( ( Ar] s 
£i{r,z,<p,t) = — - - + -- 
c s T \s 3s / 
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where the real constants to be chosen are A, B and D. 



(235) 



VII. NUMERICAL SOLUTIONS 



In this section we present numerical results. The solutions of the differential equations can 
be grouped in those which are smooth and those which exhibit some non-smooth behavior, 
such as rapid oscillations or the formation of "walls" , specially at later times. This kind of 
behavior appears when there is a lack of balance between the different terms of the equations. 
Therefore, before presenting numbers and plots, we discuss, in the next subsection, the 
conditions for finding stable solutions. 
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A. Soliton stability 



As we have seen, perturbations in fluids with different equations of state generate different 
nonlinear wave equations. Some examples are the Kadomtsev-Petviashvili (KP) equation: 

d [dF 8F d 3 F) d 2 F d 2 F , , 

^7 + + a 2^r +"3^ + a 41 ^ = 0. 236 

ox { at ox ox 6 J Oy z Oz z 

and its particular cases, such as the KdV: 

dF dF d 3 F . . 

m +a ^ + a ^ = ° (237 > 



and the breaking wave equation: 



dF dF . , 

M +a ^=° (238) 



Another example of a nonlinear wave in a dissipative system is the Burgers equation: 

dF dF d 2 F , . 

m +aiF ^ + a ^-° (239) 

where ol\ to are real constants. Having derived a particular differential equation, we 
can check whether the obtained equation is consistent with the physical picture of a small 
amplitude and long wave length perturbation propagating over large distances. We shall 
follow the analysis performed in Ref. [18]. Let us assume that the nonlinear equations listed 
above for a generic function F have a solitary wave solution with a typical large length 
L ~ l/cr (0 < <t << 1). Considering the general case, the KP equation has a dispersion 
term that is about ~ a 4 F. It must arise at a propagation distance (or equivalently 
propagation time T) D, accounted for in the equation by the term ~ a^. If both 
the dispersion and propagation terms have the same size, then T ~ l/cr 3 . Regarding the 
nonlinear term, if it has the form §~{F^) its order of magnitude is a 2 F 2 . The formation 
of the soliton requires that the nonlinear effect balances the dispersion. Hence it must have 
the same order of magnitude and F 2 a 2 = Fa 4 . Hence F ~ a 2 . We can then conclude that 
F « L « D and the above equation describes the propagation of a wave with small 
amplitude (F) and large wave length (L) which travels large distances (D). In the case 
of KP we have terms that describe the transverse evolution of the wave. We can estimate 
their sizes only if we make assumptions about the transverse length scales. In most cases the 
resulting flow is one- dimensional along the x direction with some "leakage" to the transverse 
directions. 
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B. Numerical analysis 



In what follows we apply the numerical tools developed in the appendix to several cases. 



Nuclear soliton 



We start our numerical analysis showing in Fig. [2] the solution of the linear KdV equation 

we use the analytical 



at T = 0, (G = and X = x) (204) with d 



2' 



In Fig. 



2 a 



solution (220) as initial condition. As expected this pulse propagates without dissipation 



nor dispersion: it is a soliton wave. This is the situation illustrated in Fig. [T] Any change 



in the initial condition has noticeable consequences as it can be seen in Fig. |2(b) , where 



we follow the evolution of the numerical solution of (204) for an initial pulse given by (220) 



multiplied by a factor 20. As it can be seen, the amplitude grows, the width decreases and 
secondary bumps appear propagating behind the first. 
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(a) (b) 

FIG. 2: a) The evolution of the analytic solution of the KdV equation, b) The evolution of the 
analytic solution multiplied by a factor 20. 
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In Fig. [3j we show the equivalent plot for the spherical case: Q = 1 and X = r. In 
contrast to the linear case there is a strong damping of the pulse. The dependence on the 
initial conditions is also strong. 



0.10 -, 0.7 n 




r(fm) r(fm) 



(a) (b) 
FIG. 3: Analogous to Fig. [2] for spherical coordinates. 



In Fig. 4 we show for the linear case and for the "optimal" initial condition (220) the 



evolution of the pulse with time for different temperatures. We can see that, increasing the 
temperature the pulses move faster and go farther. The same feature can be observed in the 
spherical shown in Fig. \5\ . 



Setting d — — g in the wave equations (204) and (207) we eliminate the third order 



derivative terms. The corresponding wave equations are breaking wave equations. Out of 



smooth initial perturbations, given by (220), these equations create shock waves. We can see 



this process in one dimensional Cartesian coordinates in Fig. [6] . We observe a steepening of 
the profile until the formation of the shock, followed by the dispersion of the wave. We see 
that the higher is the initial amplitude, the sooner the wave breaking and dispersion occurs. 
In Fig. [7] we fix one initial profile and study its time evolution for two different tempe- 
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FIG. 4: The panels show calculations with temperatures a) T = 20 and b) T = 70 MeV . 

ratures. 



7(a) and 7(b) show the development of a shock wave at T = 20 MeV and 70 MeV 
respectively. As it can be seen, with increasing temperatures the pulse moves faster and the 
shock formation and the subsequent dispersive breaking occurs later. For the radial case a 
similar behavior is observed. 



2. Breaking wave equation in QGP 



The KdV equation can be written as 

du du n d ( u 2 

which has the following analytical soliton solution: 



d 3 u 
dx 3 



f(x,t) 



3(u — c s 



sech 



\u - c a ) 



4/i 



(x — ut) 



In the numerical study of (208) and (209) we use the following soliton-like profile: 



Pi(x, t ) = A sech 2 



(240) 



(241) 



(242) 
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(a) (b) 
FIG. 5: Analogous to Fig. [4] for spherical coordinates 

In this equation A and B represent the amplitude and width of the initial baryon density 
pulse, respectively. 

we 



We present in Fig. pUhe numerical solution of (208) for different times. In Fig. 8 (a 



show for A = 0.075 and B = 1 fm and in 8(b)| for A = 0.35 and B = 1 fm. It is possible 
to observe the evolution of the initial gaussian-like pulse with the formation of a "wall" on 



the right side. In 8(b) the "wall" formation and dispersion occurs much earlier than in 8(a) 



due higher initial amplitude. 



In Fig. 9(a) we show the solution of (209) with the initial condition given by (242) with 
A = 0.01, B = 1 fm and T = 300 MeV. Fig. |9(b)| shows the same but with A = 0.1 and 
B = 1 fm. As in the zero temperature case, increasing the initial amplitude the breaking 



process and dispersion develops earlier. From Fig. 9(a) we can conclude that it is possible 



to find an approximate solitonic behavior even when the differential equation is not the KdV 



one. 



47 



10 20 30 40 50 

x(fm) 



5 10 15 20 

x(fm) 
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FIG. 6: a) Shock wave formation in one dimensional Cartesian coordinates, b) Initial profile 
multiplied by a factor 20. 

3. KP equation in cold QGP 



We now study the conditions in which the solution (223) must be real and therefore the 
constant h\ must be positive. Moreover, following Refs. [33j [31] we assume that a 2 + b 2 = 1 
and consider p\ a normalized perturbation [12] ; within the region (in the u — a plane) Eq. 
(223) is well defined and we can have solitons. The parameters are chosen to be: po = 1 
fm~ 3 , g = 1.15 and vtlq = 460 MeV, which imply c s ~ 0.64 [12]. The stability analysis 
can be made more rigorous with the introduction of the Sagdeev potential [331 [31] by using 



r] = £ — ut = ar + bz — d 



ut to rewrite equation (212) as an energy balance equation. 



For our present purposes the requirements in [12] are sufficient 



The plot of the soliton evolution is presented in Fig. 10 and in Fig. 11 We show a 



plot of (223) with fixed ip = 0°, a = 0.6, b = 0.8, u = 0.73 and z varying in the range 



Ofm < z < 30 fm which satisfies the soliton conditions. In Fig. 10(a) the pulse is observed 
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(a) (b) 

FIG. 7: Shock wave formation in one dimensional Cartesian coordinates for different temperatures: 
a) T = 20 and b) T = 70 MeV . 



at t = 18 fm whereas in Fig. 10(b) at t = 28 fm. From the Fig. 10(a) we can see that the 
cylindrical pulse expands outwards in the radial direction. The regions with larger z expand 
with a delay with respect to the central (z = 0) region. 



Keeping z = 1 fm fixed, we show the time evolution of (223) from t = 10 fm (Fig. 



11(a)) to t = 22 fm (Fig. |ll(b) ). The azimuthal angle varies in the range 20° < (f < 150°. 
From the parenthesis in ( |223 ) we can see that the expansion velocity grows with the angle. 
This asymmetry can be clearly seen in the figure, where the large angle "backward" region 
moves faster the small angle "forward" region. The breaking of z invariance and azimuthal 
symmetry is entangled with the soliton stability [32] and with the physical properties of the 
system (contained in the parameters hi, h,2 and c s ). 



We perform the study of the existence condition for the solution (227), which must be 



real and therefore the constant U — w must be positive. The parameters are the same: 
p = 1 fm~ 3 , g = 1.15 and m G = 460 MeV, which imply c s ~ 0.64 [32]. We also set C = 0.5 
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(a) (b) 
FIG. 8: Time evolution of a baryon density pulse for cold QGP. 

and extend the condition in Refs. [33J [M] to A 2 + B 2 + C 2 = 1 . As mentioned before, 
U > w and again, p\ is a normalized perturbation jl2] and within the region (in the U — A 



plane), (227) is well defined and we can have solitons |42j . The stability analysis can be 
performed more rigorously with the introduction of the Sagdeev potential [331 El] by using 



Ax + By + C z — Ut, to rewrite equation (213) as an energy balance equation. A simple 



example of soliton evolution is presented in Fig. 12 The plot of (227) with fixed z — 1 fm, 
A = 0.6, B = 0.62, U = 0.66 and y varying in the range Ofm < y < 50 fm. The pulse is 



observed at two times: t = 30 fm (Fig. 12(a)[ ) and t = 120 fm (Fig. |12(b) ). From the figure 
we can see that the cartesian pulse expands outwards in the x direction keeping its shape 
and form as expected. 



Analogously to the nuclear soliton, the particular case of KP (213): the KdV (214), has 
the exact soliton solution (229). In Fig. 13 we choose u = 0.8, po = 2 fm~ 3 and c 2 = 0.5. 



In Fig. 13(a), the numerical solution of (214) for (229) as initial condition is studied 
for different times. As expected, the evolution of the initial gaussian-like pulse as a well 
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(a) (b) 
FIG. 9: Time evolution of an energy density pulse at T = 300 MeV. 



defined soliton, keeping its shape and form. In Fig. 13(b), we show again the numerical 



solution of (214) for (229) multiplied by a factor 10. Now the initial soliton posses amplitude 



6.0 and starts to develop secondary peaks, which are called "radiation" in the literature. 
Further time evolution would increase the strength of these peaks until the complete loss of 
localization. 



For the other particular case (215) of the KdV (214), which we show in Fig. 14, where 



again we use (229) with the same parameters listed above. Fig. 14(a) , we have the breaking 



following by its dispersion of the initial pulse. And in Fig. 14(b), we use (229) multiplied 
by a factor 10, which we observe the anticipation of the breaking following by its dispersion 
of the initial pulse much earlier in comparison to Fig. |14(a) 
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(a) (b) 
FIG. 10: Cylindrical soliton time evolution in r — z plane. 



4- Burgers equation in hot QGP 

We apply the perturbation concept to study the radial expansion of cylindrical flux tubes 
in a hot QGP. These tubes are treated as perturbations in the energy density of the system 
which is formed in heavy ion collisions at RHIC and LHC as explained in [33]. During 
the expansion there is a "competition" between the background and the tube. In Fig. 



15(a) the QGP background expands faster and the tube is "pushed" outwards generating an 



anisotropic energy (and final particle) distribution. In Fig. 15(b) the opposite occurs: the 
tube expands faster than the background, generating a different type of anisotropy in the 
final state. In principle two and three-particle correlation measurements could distinguish 
between the two cases. Viscosity may change this picture. As shown in [13J, a strong 



viscosity could rapidly damp the tube and reduce any anisotropy. In Fig. 15 we have two 
extreme situations and something in-between may occur. With our formalism we can study 
quantitatively the evolution of the tube. 
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(a) (b) 
FIG. 11: Cylindrical soliton time evolution in r — tp plane. 



We perform the numerical analysis for equations (218) and (219) with the initial condition 
given by a gaussian pulse in i\\ 

ei(r) = Ae- r2/r ° (243) 

where the amplitude A and the approximate width r$ are parameters which depend on the 
dynamics of flux tube formation. We shall refer to r as the initial "radius" of the tube. The 
tubes are perturbations, so we expect A < 1. According to |43J and references therein, the 
transverse size of the tubes is of the order of 1 fm and in our calculations tq = 0.8 fm. We 
consider hot QGP at temperatures Tq = 150 MeV and Tq = 500 MeV treated as an ideal 



fluid {rj/s = (/s = 0) described by (219) and as a viscous fluid (rj/s = 0.16 and (/s = 0) 



described by (218) [43J. 



In Fig. 16 we show numerical solutions of (218) for a viscous fluid using (243) and we 
can observe the increasing temperature favors the tubular structure survival. 



In Fig. 17 we perform the same study for the ideal fluid (219) with (243) and we show 
that breaking with dispersion occurs. 
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(a) (b) 
FIG. 12: Cartesian soliton time evolution in x — y plane. 



By comparison between Fig. 16 and Fig. 17 we conclude that viscosity dissipates the 
breaking followed by dispersion of the pulse. The tube expands radially with a supersonic 
velocity and in less than 4 fm/c it becomes a "ring", with a hole in the middle. Moreover, 
by this time the amplitude is already reduced by a factor two and the tube (or ring) looses 
the strength to "push away" the surrounding matter |43j. 



5. KP-Burgers in hot QGP 



As an example of time evolution for the analytical solution of the cKP-B equation (217) 
we plot the time evolution of ((2321) with <p = Orad, D = 1, A = B = 0.5, T = 300 MeV 



and the viscous fluid with r]/s = 0.16 and (/s = 0. The Fig. 18 shows the time evolution of 
the analytical traveling wave. 



As it can be seen, (232) obtained depends directly of the dissipative (viscous) coefficient 



term of (217), i.e., it is a dissipation dominated solution. The Fig. 18 shows a shock wave 
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(a) (b) 

FIG. 13: a) Soliton propagation in one dimensional Cartesian coordinate and b) Initial profile 
multiplied by a factor 10 providing peaks plus "radiation" . 

propagation in a time interval of 50 fm. 



VIII. CONCLUSION 



The discovery of the quark gluon plasma in the high energy heavy ion colliders brought 
relativistic hydrodynamics to the main stage of hadron physics. Encouraged by the vigorous 
experimental program at CERN theoreticians of hydrodynamics embarked in an ambitious 
project: the calculation of observables quantities with relativistic viscous hydrodynamics. 
The measurement of two and three particle correlations may be useful to study the propa- 
gation of waves in the QGP. Among the sources of waves we have fast partons crossing the 
medium and also flux tubes formed in the initial stage of heavy ion collisions. In this work we 
have emphasized that these waves are most likely nonlinear and should be studied with the 
appropriate formalism, which, in our opinion, is the Reductive Perturbation Method. After 
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(a) (b) 
FIG. 14: a) Breaking wave and b) breaking wave anticipation. 

making a survey of relativistic hydrodynamics we have presented the RPM in a simple and 
pedagogical way. The equation of state of the two most relevant strongly interacting fluids, 
i.e., nuclear matter and the quark gluon plasma, was discussed. In both cases we have made 
an effort to give a pedagogical introduction for non experts. In both cases we have shown 
how to obtain a KdV soliton. The main responsible for the appearance of these solitons are 
higher order derivative terms in the vector fields appearing in the Lagrangian of the system. 
In QHD it is enough to relax the strict mean field approximation, in which all gradients 
vanish, and allow for spatial inhomogeneities in the vector field. A slightly more careful 
treatment of the vector field equation of motion yields the desired term, derived from the 
Laplacian V 2 V . In QCD, the situation is more complicated because the gluon is massless. 
This makes impossible a simple estimate of the corresponding quantity V 2 Aq. However a 
careful treatment of the non-vanishing vacuum condensates leads to a dynamically generated 
gluon mass, mg, which introduces a mass and a size scale and renders possible the estimate 
of the desired Laplacian and the existence of solitons in the QGP. 
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(a) (b) 

FIG. 15: The large circle represents the front view of a cylindrical portion of QGP fluid. The 
small and dark circle represents a flux tube, i.e., a cylindrical perturbation with a higher energy 
density than the background, a) The tube expands slowly and the background faster, b) The tube 
expands much faster than the background. 

Combining the equations of hydrodynamics with the equation of state and applying 
the RPM we have derived several differential equations for the perturbations in energy 
and baryon density. These equations connect properties of the waves, such as width and 
speed, with the microscopic dynamical quantities of the fluids, such as particle masses and 
couplings. Several of them have analytical solutions, which were presented and discussed. 
Some others must be solved numerically. As expected for nonlinear equations, the results 
depend very strongly on the initial conditions. An interesting finding is that, even when 
we do not have a KdV equation, in many cases the breaking wave equation has very long 
living localized solutions, which resemble to solitons. In some other cases, the initial pulses 
loose their localization and/or start to present rapid oscillations. All these features may 
manifest themselves directly or indirectly in the experimental data. The analysis made here 
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FIG. 16: Solutions for several times: a) T = 150 MeV and b) T = 500 MeV. 

is still qualitative and a closer contact with phenomenology is still to be made. For now, 
the obtained results suggest that viscosity strongly affects the propagation of perturbations 
in the quark gluon plasma. In order to confirm this statement the next step is to apply the 
RPM to the Muller-Israel- Stewart theory. 

We hope to have convinced the reader that the study of nonlinear waves in hadron physics 
is an interesting and fast moving field. This study will help to interpret and understand the 
data from the LHC. 

IX. APPENDIX: METHOD OF FINITE DIFFERENCES 

The most general form of a one-dimensional nonlinear wave equation with a second order 
dissipative term and a third order dispersive terms is given by: 




(244) 
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(a) (b) 
FIG. 17: Analogous to Fig. pi but for the ideal fluid, a) T = 150 MeV and b) T = 500 MeV. 



For a = 0, (5 and fi ^ 0, it is the usual third order nonlinear Kortweg-de Vries equation 
(KdV). For fi = 0, and a ^ , it represents the Burgers equations. Finally, when only 
(3 7^ 0, we have the breaking wave equation. In order to solve this equation numerically, 
we divide the integration region (0 < x < (n + l)h and < t < mAt) with a space step 



h and a time step At. Then, the wave function u(x,t), solution of equation (244), assume 
the discrete values w i)0 ,Wi,i, Wij-2, u i,j+i, u i,j+2, ■■■,u iin ,u i>n+ i in a given time U. 

The expansion of this solution in Taylor's series leads to: 

— ' + (s).. h+ s(£).r + a(£).* + <™ (245) 

or, alternatively, to: 

(s).. k+ a(») .." 2 -M£0.- ft3+o(ft4) (246) 

We combine these expressions conveniently to obtain finite difference expressions for first, 
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(a) (b) 
FIG. 18: Shock wave evolution: a) t = fm and b) t = 50 fm. 



second and third order centered partial space derivatives [HI 



dx 



1,3 



2h 



O x Ui, j = 2, • • • ,n - 1 



' d 2 u 
dx 2 



i j 



u i,j+l ~ ^ U i,j "I" U i,j-i 
h 2 



O xx Ui j = 2, • • • , n - 1 



and 



'03 



dx 3 



2h 3 



OxxxUi j 2, • • • , Tl 



in which we define 



i u i3 \ 



u i,j+l 



U 



\ Ui,n-l J 
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{uifl = Ui t \ = Ui t n = Uj >n +i = are the boundary values on x axe) and the operators: 



/ l/2h 
-l/2h l/2h 



o \ 









■l/2h l/2h 
-l/2h J 



( -2/h 2 1/h 2 
1/h 2 -2/h 2 1/h 2 



o ^ 





V 



o 
o 



1/h 2 -2/h 2 1/h 2 
1/h 2 -2/h 2 J 



and 



( -l//i 3 l/2/i 3 
1/h 3 -l//i 3 l/2h 3 
-l/2/i 3 l//i 3 -1/h 3 l/2h 3 



Or., 









V 



o 






-l/2h 3 1/fi 3 -1/h 3 l/2h 3 
-l/2h 3 1/h 3 -l//i 3 
-l/2h 3 1/h 3 J 



Analogously, we can expand the solution of equation (244J) around a time tf 

' du 



u '» = u '^ft). At+ l(w). (At)2+0{h3) 



to obtain 



' du 

at 



Ui+i - Uj 
At 



(251) 



(252) 



(253) 



(254) 



(255) 



When we apply (255) and the operators (247), (248) and (249) in equation (244), it can 



be represented equally in times £j or £ i+1 , i.e., the following: 

-P0 x [{l/2)Uf] - aO xx U t - fjiO^Ui 



U i+ i - Ui 
At 



or 



Ui+i - Ui 
At 



-PO x [{l/2)U 2 +1 ] - aO xx U i+1 - fiO xxx U i+l 
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(256) 



(257) 



As long as the error one makes in both situations is the same, the Crank-Nicolson scheme 
prescribes to take the mean value of these two possibilities, therefore 



U i+l =Ui- ±AtpO x [(l/2)U? + (l/2)C/f +1 ] - ^AtaO xx (U t+1 + U { ) - ^AtfiO xxx (U i+1 + U { ) 

(258) 



Given a initial condition U , U i+ i is the solution of (258) at any posterior time. However, 
this is a set of nonlinear algebraic equations. This problem becomes much more simple if 
we linearize it. We expand the nonlinear term in a Taylor's series, around the i — th time t% 
[63]: 

l„o . . fd(l/2)u 2 \ (du\ +0{Af2) 



2^ = 2^ + A \ On 



du 
dt 



which leads to 



luli = lUi+Ui{U w -U t 



(259) 
(260) 



where we have used (255). 



Replacing this result in (258) we get: 



U i+ i = Ui- ^AtpO x (UiU i+1 ) - l -AtaO xx {U i+l + UJ - ^At/iO xxx (U i+1 + U t ) (261) 



Using the matrix format of the operators O x (|251|), O xx (|252|) and O xxx (253), equation 



(261) becomes a quin-diagonal algorithm: 





1 - 2s -q itj + s -p 

qij + s 1 - 2s —qij + s —p 
p qij + s 1 -2s -qij + s -p 















in which 



p qij + s 1 - 2s -q^ + s -p 
p qij + s 1 — 2s —Qij + s 

p qij + s 1 — 2s J 



Wi+1,2 ^ 




1 r i(2 


Ui+lj-l 






u i+l,j 




r i,j 


u i+l,j+l 




r i,j+l 



\ ^i+l,n-l / 



P 



1 At 

1^ 



\ r i,n-l J 

(262) 



(263) 



1 At 

%3 = --^-^Uij^-^V 



(264) 
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and 



^1 + 2s)uij + pUij +2 - {2p + s)u i)j+l + (2p - s)u itj -i - pUij-2 (266) 



Therefore, from a given initial condition: 



U 



1 Uqo ^ 



U 0,j-1 

u 0>j 

U 0,j+1 

V «0,n-l y 



(267) 



the set of linear algebraic equations (262) can be iteratively solved, to obtain the solution 



u(x,t) of equation (244) in any posterior time ti+% is given by: 



U i+1 



1 Ui+1,2 ^ 



Ui+lj-l 

u i+l,j 
u i+l,j+l 

\ Wj+l 5 n-l J 



(26£ 



We turn now to the two dimensional extension of the KdV equation (237), the so called 



Kadomtsev-Petviashvilli (KP) equation (236): 

d f du ot\ du 2 d 3 u 1 d 2 u 
dx\dt 2 dx ° 2 dx 3 J ° 3 dy 2 

Repeating all the preceding procedure for this equation, we obtain: 



(269) 



O x U i+1 = OJJi-^Ma^U^^+^-^M^O^U^+i + U^-^Ma^OyyiU^ + U^ (270) 
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in which we define, 



Ui 



1 u i)2 ,k ^ 



U i,j+l,k 



k = 1,2,...,/ 



(271) 



i Ui, n -l,k J 

{u>i,j,o = u i,j,i+i — are the boundary values in the y direction) and the operators [61] : 

( 6/h 4 -A/h 4 l/h 4 ■■• \ 
-A/h 4 6/h 4 -A/h 4 l/h 4 ■■• 
l/h 4 -A/h 4 6/h 4 -A/h 4 l/h 4 ■■■ 



O, 








l/h 4 -A/h 4 6/h 4 -A/h 4 l/h 4 
l/h 4 -A/h 4 6/h 4 -A/h 4 
1/h 4 -A/h 4 6/h 4 j 



(272) 



is the finite differences fourth order centered partial space derivative operator, and 

( -2/h* l/h 2 y ••• o \ 
l/h\ -2/hl 1/hl ■■■ 



o 



yy 



(273) 



o • • • i/hi - 2 /K l / h l 

\ o ••• o i/h 2 y - 2 IK) 

in which h y is the space step in the y direction. 

Replacing the matrix representation of the operators in equation (270) we get: 

/ 



c i,j,k di,j,k a 

&i,j,k CiJ,k di,j,k & 
O bi,j,k Ci,j,k dijk O 














^ ™i,j,k ^i,j,k dijfc d 
Ct bi,j,k ^i,j,k dijk 

0a b i>jtk Ci >j>k J 



( 



\ 



u i+l,j-l,k 

u i+l,j,k 
u i+X,j+l,k 

\ Ui+l,n-l,k } 



( 



\ 



e i,2,k 

@i,j—l,k 
&i,j,k 

\ &i,n—l,k j 



k = 1,2,...,/ 



(274) 
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in which 

At 

a = (275) 

hj,k = 1 - 4a + ^—^u i:j+li k (276) 

2a 3 hAt 2 ai At 
Ci,j,k = 6a — — m^k (277) 

ot\At 

dij,k = -1 - 4a H — Wij-i.fe (278) 

2«3aAt 

ejj,fc = -au i)j+2 ,k + (1 + 4a)w iii+ i ifc + ( — 7^ 6a)Uij,k + (-1 + 4a)M ij _i ifc - au i; j- 2> k 

2a 3 hAt 

(Wj+i,i,fc+i + Wi + ij 5 fc_i + Uij^k+i + Uij^k-i) 



hi 



The algorithm (274) represents / sets of linear algebraic equations. As long as each one 
of these sets are self-consistent, they must be solved iteratively, until the answer converges 
to the solution (U l+1 )f after / iterations. To stop the iterations, we can use the criterion 

11(^1)/ -(^)/-i 

ll(^+l)/-lll <e ^ Z(J ) 

In what concerns the stability of this numerical method, the conservation of some quan- 
tities such as: 

rco rco 

Pit) = / / udxdy (280) 



— OO J —OO 



j poo poo 

E[t) = - / u 2 dxdy (281) 



and 

E(t) = 

A J —CO J —CO 

are frequently used as a criterion to verify its reliability. In refs. [65] and [66], this criterion 
is used to show analytically that these numerical methods based on the Crank-Nicolson 
scheme are unconditionally stable. 
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